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Tlie divergence of tlie tliermal conductivity in the thermodynamic limit is thoroughly inves- 
tigated. The divergence law is consistently determined with two different numerical approaches 
based on equilibrium and non-equilibrium simulations. A possible explanation in the frame- 
work of linear-response theory is also presented, which traces back the physical origin of this 
anomaly to the slow diffusion of the energy of long-wavelength Fourier modes. Finally, the 
results of dynamical simulations are compared with the predictions of mode-coupling theory. 
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An extremely idealized, but physically meaningful, model of an insulating solid is a set of iV atoms of mass 
m, arranged on a 1-d lattice with spacing a and coupled by nonlinear forces. Denoting with qi the displacement 
of the Z-th particle from its equilibrium position la, the corresponding Hamiltonian reads 

(1) 

where pi — mqi . If the chain is put in contact at its boundaries with two heat baths at different temperatures 
Tq and Tq -I- AT, a nonequilibrium stationary state arises, characterized by a non vanishing average heat flux 
J. The microscopic definition of J is 

I I 

where fi = —V'{qi — qi-i) is the interaction force and ji represents the local flux at site I (i.e. the sum of the 
fluxes of potential energy from its neighbours). The lattice thermal conductivity k can be defined through the 
Fourier law, 

W = -4 , (3) 

where x — la is the coordinate along the chain and (•) denotes a time average in the stationary regime. 
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The task of statistical mechanics is to explain such a macroscopic law, starting from the microscopic interac- 
tions. Long time ago, Peierls proposed a theory based on the phonon scattering mechanism in the perturbative 
and close-to-equilibrium limits B. The common belief was thereby that a microscopic theory able to include 
the effects of disorder and/or nonlinearity should ensure the validity of Eq. (|^) also in the case of strong 
anharmonicity and for systems far from equilibrium. 

Nevertheless, the rigorous treatment of isotopically disordered harmonic chains showed that disorder alone 
is not sufficient for obtaining a reliable model: the thermal conductivity diverges as k (x 

iVi/2 g. This result 

enforced the necessity to verify whether nonlinearities alone can eliminate such an anomaly. Recently, the 
problem has been reconsidered by resorting to nonequilibrium molecular-dynamics simulations. Among the 
many attempts, the most convincing results in favour of a finite k were obtained only for rather artificial 
systems, namely the ding-a-ling models [|||^. 

For more realistic systems, like the diatomic Toda chain a finite conductivity was found for small enough 
mass ratio while the analysis of model (|l|) with the Fermi-Pasta-Ulam (FPU) potential 

v{y) = + \9y^ , (4) 

gives evidence of a power-law divergence k oc for ^ 1 . We have verified that this result is robust 
against changes in the number of thermostatted atoms, in the boundary conditions and in the type of heat baths 
(time-reversible as well as stochastic ones) . More precisely, simulations reported in Q| with a large temperature 
gradient (AT = 128 and Tq 24) reveal a crossover from /3 « 0.45 for N < 250 to /3 « 0.38 for larger N, while 
new simulations made with AT = 16 and Tq = 80 (Fig. 1) yield (3 = 0.37 up to TV = 2048. 

For sufficiently small AT and large N, the system is close to an equilibrium state at temperature T = To + ^. 
One can then rely on the Green-Kubo formula [0| 
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Cj{t)dt (5) 



where Cj{t) = N{J{t)J{0)) is the correlation function at equilibrium. We have numerically computed Cj{t) for 
the FPU chain with Born- Von Karman boundary conditions, and initial conditions sampled according to the 
canonical Gibbs measure. In Fig. 2, we report the power spectrum S'(w) of the heat flux J for N = 512. In the 
low-frequency region there is a power-law divergence as w"" '^^, which corresponds to a decay of Cj(t) as t~°-^'^. 

A crucial information for comparing these results with those obtained with heat baths comes from the spatio- 
temporal correlation function Cj{l, t) = (j;(t)jo(0)) of the local flux. Indeed, the peaks shown in Fig. 3 indicate 
that energy propagates at constant velocity Vp. Therefore, one can estimate the dependence of k on A^ from the 
asymptotic behaviour of Cj{t), by restricting the integral in (||) to times smaller than the "transit time" N/vp. 
This amounts to ignoring all the contributions from sites at a distance larger than A^. With the above estimate 
of Cj, one obtains that k cx A^*^'^^, in agreement with the direct estimate obtained from the non-equilibrium 
simulations in Fig. 1. 

Despite the remarkable consistency of the numerical results, the scenario is quite unsatisfactory on a theoreti- 
cal ground, especially because there is neither a clue about the very reason of the observed anomalous transport, 
nor an explicit estimate of the measured exponent. In the following, we present a possible explanation based 
on a kind of hydrodynamic approach with long-wavelength modes interacting with a stochastic-like reservoir. 

In fact, it has been observed that even at high energies, where the FPU-system is definitely nonintegrable, 
the energy of its long- wavelength Fourier modes diffuses very slowly (with respect to their period). Accordingly, 
they may effectively act as undamped transport channels, a feature that reflects itself in the relaxation of 
fluctuations as measured by Cj{t). Our goal is now to obtain an estimate of the asymptotic behaviour of Cj{t), 
by describing the relaxation close to equilibrium of the Fourier- mode amplitudes Ak and Bk, defined by 



N/2 



^ , 2TTka\ ^ . / 27rA;a, 
Ak cos I — I + Bk sm I — 



N \ N 



(6) 



where we have made the assumption that the zero-mode Aq is identically zero. The customary approach would 
be to write stochastic equations for Ak and Bk- The general strategy involves projection on the subspace of such 
(slow) variables and yields linear Langevin equations with memory terms H. If we assume that a separation of 
proper time scales is possible, the former reduce to Markovian equations |[ll[| , 

Ak + ikAk + CjlAk = a ; Bk+ ikBk + CjlBk = rik (7) 
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where 7^ are phenomenological relaxation rates, ^/c and rjk are mutually independent, Gaussian, random and 



white processes |12 . Nonlinear effects are taken into account by renormalization of the bare dispersion relation 
as 

LOk = a u)k = a ■ 2cjo sin f j , (8) 

i.e. by renormalizing the sound velocity from v — acoo to v — av, where, for the parameters of Fig. 1, 
a = v/[l + 2(1 + 0.72 e)V2]/3 and e{T) is the internal energy 0. 

A first confirmation that transport is basically determined by the low fc-modes comes from the observation 
that V is remarkably close to Vp: for the parameter values considered in Fig. 3, v = 2.38, while Vp = 2.47. 

Let us now split V into its harmonic and anharmonic parts and consequently write the flux (H) as J = Jh + Ja- 
From Eq. (^) one obtains 

N/2 

Jh = ^^CkOJk(^AkBk- AkBk') , (9) 

k=l 

where Ck = duok / dk is the bare phase velocity. For a strongly anharmonic system, like the one we consider here, 
we do not expect J a to be negligible. Nevertheless, we shall now argue that the leading asymptotic behaviour 
of Cj can be estimated by computing the autocorrelation of Jh alone. We limit ourselves to the FPU case (0), 
for which J a is a sum of fourth order terms coupling modes with indices k,k' ,k" , k'" . By virtue of Eq. (|^), each 
low-fc component oscillates approximatively at a frequency that is an algebraic sum of Wfc, Wfc', Wfc" and ujk'"- 
It is then reasonable to assume that rapidly oscillating terms will be relevant only for short times, while the 
asymptotic dominant contributions originate from the resonant terms. This fact, together with selection rules, 
require that only terms fulfilling the constraints 

k + k' + k" + k"' ^ (10) 



must be retained. For small fc, we can linearize the dispersion relation in the second of Eqs. (10) obtaining 



■'--'"s^j^E'-n^+B^) . (II) 

k 

In order to evaluate the contribution of J a to Cj, let us first notice that the former is the product of two 
time-dependent terms: Jh and the sum S ~ S'^K^l + ^fc)- given fc-contribution in the expression of 

Jh (see Eq. (^)), is correlated with only one (with the same k) out of the N addenda of the sum in Eq. ([ll|). 
Accordingly, in the thermodynamic limit, we can consider S as uncorrelated with Jh- Moreover, since S is the 
sum of uncorrelated positive terms, we can substitute it with its average value and write 

J ^ Jh(i + -4^U{T)] , (12) 



where we have also taken into account energy equipartition in its generalized form, i.e. mw^(A^) — muj'^{Bf^) = 
U(T), for some function U independent of k. We numerically checked the reliability of our approximations by 
first comparing expression (^l|) with Ja and then by verifying that the autocorrelation functions of Jh and J 
are proportional. 

Let us compute the autocorrelation of Jh- Using Eqs. (|^), one finds that the new variables Wk = AkBk—AkBk 
satisfy the Langevin equation 

Wk = -jkWk + Ck ■ (13) 
In the limit of small jk, Ck is a Gaussian and delta-correlated random process, so that for large N 

{JHit)JHm ^ ^kBTUiT)J2ck^e-''^' ^ ^-ksTUiT) T^" dfc c2(fc) e^^^* . (14) 
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Due to the conservation of total momentum, one expects 7^ to vanish for vanishing k. By drawing an analogy 
with the known nonanalytic behaviour of dispersion relations of low-dimensional fluids p3[ |, we may assume 
7(fc) « vk^^, to leading order in k [Q. We also suppose that /i > 1 consistently with the observed fact that 
dissipation occurs on time scales much longer than the proper period. Since the conduction properties are 
determined by low-fc modes, we can approximate Eq.(^) by expanding c{k) around fc = so that 



The power-law decay of Cj{t) implies that k should diverge as N^^^^f^. 

The key point is then the evaluation of fi. In the case of fluids, this is generally accomplished by resorting 
to a self-consistent mode-coupling approximation jl^ . In the present context, it is reasonable to conjecture, on 
the basis of Eqs. (|^), that the behaviour of low-fc modes can be assimilated to that of hydrodynamic modes in 
a Id fluid, at least at such high temperatures (see pT[ for a more detailed analytical discussion). One should 
then expect /i = 5/3 thus yielding Cj{t) oc and k oc N"^/^. This prediction is compared with the 

numerical results in the inset of Fig. 2. The discrepancy with the best-fit value (0.37) can be definitely attributed 
to finite-size effects. As long as the above conjecture proves correct, our numerical results provide the most 
striking confirmation of the mode coupling predictions. 

The generality of our arguments suggests that any nonlinear 1-d model with an acoustic spectral component 
should exhibit the same anomaly. It has however to be admitted that the mentioned Toda model in which an 
optical branch is also present is qualitatively different. 

Finally, we observe that there is no reason for expecting the same phenomenon not to occur also in higher 
dimensions, independently in each acoustic branch. On the other hand, the exponent /i may well depend on 
d and, moreover, the integral in Eq. (^4|) has to be weighted with the phononic density of states, which, for 
small k, vanishes as k"^^^ . One can conjecture that such a mechanism does not introduce any divergence in the 
conductivity for d = 3, while a logarithmic divergence should be observed for d = 2. This prediction provides a 
unified explanation for both the anomalous transport of 1-d systems and the finite conductivity of 3-d lattices. 
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FIG. 1. Divergence of the thermal conductivity with the system size for the FPU chain (m = uio = 1, 
g = 0.1) as obtained from nonequihbrium molecular dynamics (see Refs. [^,^ for details). Statistical errors 
are displayed only when larger than the symbols' size. 




FIG. 2. Power spectrum S{u)) (in arbitrary units) of the global flux for the FPU chain (same parameters 
as in Fig. 1) with A*' = 1024 at temperature T = 110.7. The curve results from an average over 1400 
independent initial conditions. A blow-up of the low- frequency region is reported in the inset: the straight 
line corresponds to the l/oJ^^^ behaviour. 
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FIG. 3. The (normalized) correlation Cj {I, i) of the local flux for the FPU model with N = 1024 (same 
parameters as Fig. 1). 
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